Dynamical Instability of a Rotating Dipolar Bose-Einstein Condensate. 
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We calculate the hydrodynamic solutions for a dilute Bose-Einstein condensate with long-range 
dipolar interactions in a rotating, elliptical harmonic trap, and analyse their dynamical stability. 
The static solutions and their regimes of instability vary non-trivially on the strength of the dipolar 
interactions. We comprehensively map out this behaviour, and in particular examine the experi- 
mental routes towards unstable dynamics, which, in analogy to conventional condensates, may lead 
to vortex lattice formation. Furthermore, we analyse the centre of mass and breathing modes of a 
rotating dipolar condensate. 
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In recent years a considerable amount of exp erimental unit vector e is given by, in the notation of Ref. [T| 



in recent years a considerable amount ot exp erimental 
0, H and theoretical % 0, S, 0, i, @, El ED] work 
has been carried out on dilute Bose-Einstein condensates 
(BECs) in rotating anisotropic traps. Where short-range 
interactions dominate, a vortex lattice forms when the ro- 
tational frequency (ft) of the system is w 0.7cu± (where 
luj_ is the trapping frequency perpendicular to the axis of 
rotation). Insight into the mechanism of vortex forma- 
tion can be gained by noting that 0.7lu± closely coincides 
with the frequency at which certain hydrodynamic sur- 
face excitations become unstable [4, 5]. Through compar- 
ison with experimental results [l|, |2], |6[ and numerical so- 
lutions of the Gross-Pitaevskii equation (GPE) [f| [13, [ll| 
such instability has been directly related to the formation 
of a vortex lattice. 

The above results apply to conventional BECs com- 
posed of atoms of mass m with short-range s-wave inter- 
actions, parameterised via g = ^Kn 2 a/m, where a is the 
s-wave scattering length. However, a recent experiment 
has formed a BEC of chromium atoms with dipolar in- 
teractions [12(| . This opens the door to experimentally 
study the effect of dipolar interactions in BECs. Parallel 
theoretical work, using a modified GPE, has studied the 
effect of such long-range interactions on the ground state 
vortex lattice solutions [13, [13, EH • However, the route 
to generating such states has not been explored. For this 
purpose we solve the hydrodynamic equations of motion 
for a dipolar BEC in rotating anisotropic harmonic traps. 
We show that the solutions depend on both the strength 
of the dipolar interactions, e e id, and the aspect ratio of 
the trap, 7 = lu z /luj_, in stark contrast to conventional 
BECs where they are independent of both the strength 
of the interactions and 7 [3, 0] ■ In addition we evaluate 
the dynamical stability of our solutions, showing that the 
region of £1 for which the solutions are stable can be con- 
trolled via both Edd and 7. By analogy to conventional 
BECs [3, EL HI , one may expect these instabilites to result 
in vortex lattice formation in dipolar BECs. 

Consider a BEC with long-range dipole-dipole interac- 
tions. The potential between dipoles, separated by r and 
aligned by an external electric or magnetic field along a 
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For low energy scattering of two atoms with dipoles in- 
duced by a static electric fiel d E = Ee, the coupling 
constant C d d — E 2 a 2 /e [l7l Il8j|. Alternatively, if the 
atoms have permanent magnetic dipoles, d m , aligned in 
an external magnetic field B_ = Be, one has Cdd = Mo^m 
[l9j |. Denoting p as the condensate density, the dipolar 
interactions give rise to a mean-field potential 



®dd(r) = / d 3 r'U dd {r-r!)p{r') 



(2) 



which can be included in a generalized GPE [l8|, 
for the BEC. In the Thomas-Fermi (TF) regime the 
GPE describing a static dipolar BEC in a harmonic trap- 
ping potential V(r) — m,(oj 2 x x 2 + uj 2 y 2 + cu 2 z 2 )/2 is 

1" = y K^ 2 + u 2 y y 2 + u 2 z z 2 ) + gp(r) + $ dd (r), (3) 

where /1 is the chemical potential. For ease of calculation 
the dipolar potential <&dd(L) can be expressed in terms of 
a fictitious 'electrostatic' potential <j>(r) [HJ 
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(r) + 5 fp(r) 
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where <j>(r) = J a? r' p{r^) / (Air \r — r'\) and e dd = 
Cdd/Sg parameterizes the relative strength of the dipolar 
and s-wave interactions. Self-consistent solutions of Eq. 
(|3|) for p(r), <p(r) and hence $dd(r) can be found for any 
general parabolic trap, see Appendix A of Ref. [13] • 

We consider atoms trapped in a harmonic potential 
rotating at a frequency f2 about the z-axis. In the mean 
field approximation the evolution of the condensate field, 
ip(r, t), is described by the time-dependent GPE. Writing 
the condensate field in terms of a density p(r) and a phase 
S(r) and neglecting the quantum pressure we obtain the 
conventional superfluid hydrodynamic equations 



dp 
dt 



+ V-[p(»-Oxr)] = 
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l+vf^ + ^ + ^-tnxrjUo, (6) 
at \ I m m m J 

where v — (h~/m)VS is the fluid velocity field, in the 
laboratory frame, expressed in terms of the coordinates in 
the rotating frame. The stationary solutions of Eqs. ([5]) 
and © are obtained by imposing the conditions dp/dt — 
and dv/dt — 0. We look for solutions of the form 0, 3 
v = a(yi + xj), where a is to be determined. We can 
combine this with Eq. © to obtain 

» = y + tfv 2 + ^Iz 2 ) + gp(r) + <P dd (r) (7) 

where U)\ = uj 2 + a 2 — 2afi and uj 2 = uj 2 + a 2 + 2afi. 
The form of Eq. ([7]) is identical to Eq. (O . Hence we can 
use the methodology presented in Ref. [22[ to calculate 
&dd{L)- An exac t solution of Eq. ([7]) is given by 
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where uq is the central density which is given by nor- 
malization to be no = 15 AT/ (8nR x R y R z ). Following the 
results presented in Appendix A of Ref. [22j the dipole 
potential for a polarizing field aligned along the z-axis is 
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da 



o (1 + a) (nt + a) ^[k\ + a) («* + a) (1 + a) 
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with k — x,y, z, Kk ~ and 



0o 



da 



/o (1 + *) yj + a) (nl + a)(l + a) 

Thus we can rearrange Eq. (J7J) to obtain the density 

_ fi— 3geddnoK x K y f3 - y {Co x x 2 + w a y 2 + uj z z 2 ) 
9 g(l-e dd ) 

3e dd K x K y f3 x { y} uj 2 /(2(), uj z 



(11) 
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where H>. 
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UJ 2 , , 



^0z 



uj 2 (1 - 9£ddK x K v 0z) /(2C) and C = l-e^d 
Comparing the a; 2 , y 2 and z 2 terms in Eq. ((5|) and 
Eq. (JTSJ) we find the three self-consistency relations: 
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and R 2 = ^jtC- Now using Eq. $5$ in conjunction with 
Eq. (fT2"]) we find the following expression for the station- 
ary solutions to Eq. 

= (a + fi) [u 2 x - \e d d^f^0, 



+ (a - fi) \u>, 



3 U 2 x K x K yl 2 
-£dd ^ P; 



(14) 




FIG. 1: (a) The irrotational fluid velocity, a, as a function of 
the trap rotational frequency, fl, as obtained from Eq. (I14|l . 
for 7 = 1, e = and edd = (solid curve), e^d = 0.25 (short 
dashed curve), Edd = 0.5 (long dashed curve), Edd = 0.75 
(dash dotted curve) and Edd = 0.99 (dotted curve), (b) The 
bifurcation point fit versus 7 for different dipolar interactions 
strengths; Edd increases, in the direction of the arrow, from 
to 0.9 in steps of 0.1, for the first ten curves, the lowest curve 
is for edd = 0.99. 

In the limit e dd — the solutions of Eq. l[i"4"|) are inde- 
pendent of g and 7. However, for Edd 7^ the solutions 
to Eq. (fl"4"|) are dependent on both the strength of the 
dipolar interactions and the aspect ratio of the trap. 

Introducing the parameter e = (lu 2 — uj 2 )/{lj 2 +uj 2 ) to 
define the anisotropy of the trap, we evaluate Eqs. (fT3"|) 
and (fT4]) self-consistently to determine the static so- 
lutions of the hydrodynamical equations of motion in 
the rotating frame. Figure 1(a) shows the solutions to 
Eq. (fT4"|) for various values of Edd with 7 = 1 and e = 0. 
For Edd = (solid curve) we regain the results of Refs. 
0, 0| with a bifurcation point at fib = uj x /\/2 which 
exactly coincides with the vanishing of the energy of the 
quadrupolar mode in the rotating frame. For fi < fi;,, one 
solution, corresponding to a = 0, is found. For fi > fi;,, 
three solutions appear, a = and a = ±y / 2Q 2 — lu 2 /uj x 
3 . The two additional solutions are a consequence of the 
quadrupole mode being excited for fi > ui x /\2. Actu- 
ally, it is a remarkable feature of the pure s-wave case that 
these solutions do not depend upon g. This is because in 
the TF limit surface excitations with angular momentum 
M = hqiR, where R is the TF radius and qi is the quan- 
tized wave number, obey the classical dispersion relation 
ujf = (qi/m^nV involving the local harmonic poten- 
tial V = mu! 2 R 2 /2 evaluated at R [23j. Consequently 
u>i = VluJ x , which is independent of g. However, in the 
case of long-range dipolar interactions the potential &dd 
of Eq. (|4]) gives non-local contributions, breaking the sim- 
ple dependence of the force — W upon R 16]. Thus, we 
expect the resonant condition for exciting the quadrupo- 
lar mode, i.e. fi& = lui/1 (with I = 2), to change with 
Edd- In Fig. 1(a) we see that this is the case: as dipole 
interactions are introduced, our solutions change and the 
bifurcation point (fib) moves to lower frequencies. 

In contrast to the s-wave case, not only the magni- 
tude of the dipolar coupling Edd, but also the shape of 
the BEC determines the potential <&dd- For an oblate 
( K x,y > 1) BEC, more dipoles lie side-by-side, thus giv- 
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ing a net repulsive interaction, in comparison to the pro- 
late (n x ,y < 1) case where a majority sit end-to-end, in 
which configuration the dipolar interaction is attractive. 



In the extreme limits of k x . 



and k x . 



the 



angular dependence of the interactions plays no role and 
the gas behaves conventionally, but in the intermediate 
regime the role of k XiV and hence the aspect ratio of the 
trap is important. In Fig. 1(b) we plot Q}, as a function 
of 7 for various values of Edd- For Edd = we find that 
the bifurcation point remains unaltered at f2(, = uj x j \[2 
as 7 = lj z /uj x is changed 0,13]. As Edd is increased the 
value of 7 for which Jlj, is a minimum changes from a 
trap shape which is oblate (7 > 1) to prolate (7 < 1). 

Consider now the effect of finite trap anisotropies 
(e > 0). In Fig. 2(a) we have plotted the solutions 
to Eq. (fnj for various values of Edd with 7 = 1 and 
e = 0.02. As in the case without dipolar interactions 
0, |3] the solution a = is no longer a solution for all f2. 
The effect of introducing the anisotropy, in the absence 
of dipolar interactions, is to increase the bifurcation fre- 
quency fib- Turning on the dipolar interactions, as in the 
case of e = 0, reduces the bifurcation frequency. 

We now analyze two procedures for generating an in- 
stability by tracing different paths on Fig. 2(a). 
Procedure I: fl is fixed at $1 > ilb(e = 0) and the trap 
anisotropy is adiabatically turned on. Following analy- 
ses for conventional BECs 0, [H, 0] we find that as e is 
increased adiabatically, from zero, the a = solution 
moves to negative values of a and the BEC follows this 
route. However, as e is increased further the edge of the 
lower branch fi&(e) shifts to higher frequencies. At some 
critical value of e, f2&(e) = f2, the lower branch ceases 
to be a solution. In this manner the evolution of the 
static solutions induces instability, which in conventional 
BECs has been experimentally and theoretically linked 
to vortex lattice formation 0,0]. As the dipole interac- 
tions are increased the bifurcation frequency is reduced 
and the range of SI for which this type of instability can 
occur increases from [uj x /\/2,uj x ] to [0.5^, uj x \. In addi- 
tion, dipolar interactions increase the value of e for which 
lower branch solutions exist. 

Procedure II: e is fixed and SI is introduced adiabatically. 
In this case the BEC will follow the upper branch of the 
solutions. For conventional BECs, a vortex lattice will 
form 0, when the upper branch solutions (a > 0) to 
Eq. (TTJJ become dynamically unstable. Below we gener- 
alize the analysis of Ref. [4| to examine the dynamical 
stability of the solutions to Eq. (TH|) . 

Consider small perturbations in the BEC density and 
phase of the form p = pq + Sp and S = Sq + SS then, 
via Eqs. ([5J O the dynamics of such perturbations can 
be described, to first order, as 
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where K = 
v — Q x r. 



-3^ / dxdydz / (iir \r' -r\) - 1 and w c = 
As in Ref. we consider a polynomial 



FIG. 2: (a) a vs. ft, as in Fig 1(a) but for e = 0.02. (b) The 
positive imaginary eigenvalues of Eq. (15) as a function of Edd 
for Q — 0, e = 0.02, 7 = 1 and n = 2. The three solid curves 
[two at Im(A) = LQ X and one at Im(A) = lo x (1 + €) ' 5 = u v ] are 
the frequencies of the center of mass modes of the BEC. The 
dashed curves are the frequencies of the breathing modes of 
the BEC. (c) The positive imaginary eigenvalues of Eq. (15), 
as a function of S7 for Edd = 0.5, e = 0.02, 7 = 1 and n — 2. 
The solid (dashed) curves are the frequencies of the center of 
mass (breathing) modes of the BEC. (d) The maximum posi- 
tive real eigenvalues of Eq. (15) (solid curves), as a function of 
ft, for e = 0.02, 7 = 1, n = 3 and E d d = 0, 0.2, 0.4, 0.6, 0.8, 0.95 
and 0.98; Edd increases in the direction of the arrow. The short 
and long dashed curves are additional positive eigenvalue so- 
lutions for Edd ~ 0.95 and 0.98 respectively. 

ansatz, of order n in the coordinates x,y,z and evalu- 
ate the evolution operator for the perturbations. If one 
or more of the eigenvalues, A, has a positive real com- 
ponent the stationary solution is dynamically unstable. 
However, imaginary eigenvalues correspond to stable os- 
cillatory modes of the system [23] ■ Below we consider 
both the stable and unstable modes of the upper branch 
static solutions, for 7=1 and e = 0.02. 

In Fig. 2(b) we plot the positive imaginary eigenval- 
ues of Eq. (15) for n — 2 as a function of the dipolar 
interaction strength with $1 = 0. As expected we find 
three modes associated with center of mass oscillations 
[solid curves in Fig. 2(b)], two at Im(A) = w x = w z 
and one at Im(A) = w y . These modes are independent 
of the strength of the dipolar interactions since they do 
not alter the shape of the BEC. The higher frequency 
modes (dashed curves) are associated with the breathing 
modes of the system [2J] . Since these modes do alter the 
shape of the BEC and the resulting dipolar interactions 
we find that they are dependent upon Edd- However, the 
breathing mode at Im(A) — uj x ^J1> is associated with a 
perturbation in x, y and z which is equivalent to a uni- 
form re-scaling of the density and as such the frequency 
of this mode is almost independent of the dipolar interac- 
tion strength, as can be seen in Fig. 2(b). Fig. 2(c) shows 
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how these modes shift as a function of 0, for Edd = 0.5. 
Under the rotation the motion in x and y is coupled. 
Thus the frequency of the center of mass modes (solid 
curves) is shifted in the x — y plane (23j whilst remaining 
constant, at u> z , in the z direction. Again these modes are 
independent of the strength of the dipolar interactions. 
Conversely, the breathing modes (dashed curves) have a 
relatively complex dependence on Sdd and the rotation 
frequency. 

Finally we consider the real positive eigenvalues of Eq. 
(15), associated with regions of instability for the up- 
per branch static solutions. In the limit of Edd = we 
reproduce Fig. 2 of Ref. [1], with the solutions being 
unstable in the range [0.78^, u x ] for e = 0.02. In Fig. 
2(d) we have plotted the real positive eigenvalues, Re(A), 
of Eq. (15), as a function of SI for various values of Sdd 
with n = 3. For higher values of Sdd [0.95 and 0.98 in 
Fig. 2(d)] there can be more than one real positive eigen- 
value, thus we define the region of instability as the range 
over which max Re(A) > 0], as shown by the solid curves 
in Fig. 2(d) [25 1. As the dipolar interaction strength is 
increased the lower bound in SI for the unstable region is 
reduced. For example, for Edd — 0.6 the range of rotation 
frequencies where the upper branch solution is unstable is 
[0.75ui x ,iv x ], this increases to [0.67cu x , w x ] for Edd = 0.98. 

By calculating the static hydrodynamic solutions of a 
rotating dipolar BEC and studying their dynamical sta- 
bility, we have predicted the regimes of instability of the 



condensate. In general we find that the bifurcation fre- 
quency, Qb, decreases in the presence of dipolar interac- 
tions. Thus, for a fixed SI [at SI > ilb(e — 0)] and an 
adiabatic increase in e, the critical anisotropy at which 
we expect instability to occur will be higher than for a 
conventional BEC. Furthermore, we find that the size of 
this shift not only depends on the strength of the dipo- 
lar interactions but also on the aspect ratio of the trap, 
with the maximal shift being for 7 < 1 (edd — * 1)- For a 
fixed anisotropy and an adiabatic increase in Q we find 
that as Sdd is increased the lower bound on the rotation 
frequency at which a rotating dipolar gas will be unsta- 
ble to perturbations is decreased. In conventional BECs 
these instabilities have been related to vortex lattice for- 
mation [j|. This occurs, primarily, because the instability 
disrupts the BEC at an SI which is greater than the ro- 
tation frequency at which it is energetically favorable to 
have a vortex state |26j. However, in a prolate trap the 
rotational frequency at which it is energetically favorable 
to form a vortex in a dipolar BEC grows rapidly as Sdd 
is increased [27| and can exceed the frequency at which 
we expect an instability to occur. The final state under 
these circumstances warrants further investigation. 
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